from IPython import get_ipython
Famoso conjunto de imagens de dígidos manuscritos, dividido em $60000$ imagens de treinamento e $10000$ imagens de teste. Todas as imagens são em escala de cinza com tamanho $28 \times 28$. Existem 10 classes, correspondentes aos dígitos de 0 a 9.
Este dataset está disponível em diversos lugares. Abaixo, vamos usar a versão disponível na biblioteca Keras [1]. Informações adicionais estão disponíveis no site oficial.
from keras.datasets import mnist
(x_train, y_train), (x_test, y_test) = mnist.load_data()
print(x_train.shape, y_train.shape)
print(x_test.shape, y_test.shape)
from IPython import get_ipython
Famoso conjunto de imagens de dígidos manuscritos, dividido em $60000$ imagens de treinamento e $10000$ imagens de teste. Todas as imagens são em escala de cinza com tamanho $28 \times 28$. Existem 10 classes, correspondentes aos dígitos de 0 a 9.
Este dataset está disponível em diversos lugares. Abaixo, vamos usar a versão disponível na biblioteca Keras [1]. Informações adicionais estão disponíveis no site oficial.
from keras.datasets import mnist
(x_train, y_train), (x_test, y_test) = mnist.load_data()
print(x_train.shape, y_train.shape)
print(x_test.shape, y_test.shape)
Vamos visualizar algumas das imagens
import matplotlib.pyplot as plt
get_ipython().magic(u'matplotlib inline')
fig, ax = plt.subplots(2, 3, figsize = (18, 12))
for i in range(6):
ax[i//3, i%3].imshow(x_train[i], cmap='gray')
ax[i//3, i%3].axis('off')
print('label:', y_train[i])
plt.show()
Por enquanto, não iremos aplicar nenhuma técnica de aprendizado diretamente sobre as imagens (pixels). Vamos utilizar alguns atributos que iremos extrair a seguir.
No livro Learning from Data [2], um dos atributos que os autores extraem como exemplo do MNIST é a Simetria horizontal.
Seja a assimetria definida como a diferença absoluta média entre os pixels da imagem original e da imagem refletida horizontalmente, a simetria é simplesmente o inverso da assimetria.
import numpy as np
def simetry(image):
# A operação abaixo inverte a ordem das colunas da imagem
reflected_image = image[:, ::-1]
return -np.mean(np.abs(image - reflected_image))
Outro atributo extraído como exemplo no livro Learning from Data é a intensidade média dos pixels. Este atributo está associado a proporção da imagem ocupada pelo dígito. Por exemplo, a intensidade de uma imagem com o dígito $1$ é menor do que a de uma imagem com o dígito $2$ out $5$.
import numpy as np
def average_intensity(image):
return np.mean(image)
Agora, aplicamos as funções descritas anteriormente sobre as imagens do MNIST. Note que, na forma de imagens, os dados eram representados por $28 \times 28 = 784$ atributos, enquanto que agora, eles são representados por apenas $2$ atributos.
É de se esperar que um certo grau de informação se perca neste mapeamento. Em situações reais, escolhemos as features a serem extraídas de modo a preservar o máximo possível de informações úteis.
import numpy as np
# A função abaixo converte uma imagem em uma lista de atributos,
# usando as funções definidas acima
def convert2features(image):
return np.array([average_intensity(image), simetry(image)])
# Aplica a conversão a todas as entradas do dataset
x_train_features = np.array([convert2features(image) for image in x_train])
x_test_features = np.array([convert2features(image) for image in x_test])
# Ajusta a escala das features. Utilizar multiplas features com escalas diferentes pode ser problemático
for i in range(x_train_features.shape[1]):
avg = np.mean(x_train_features[:, i])
stddev = np.std(x_train_features[:, i])
x_train_features[:, i] = (x_train_features[:, i] - avg) / stddev
# (Sim, as features no conjunto de teste são ajustadas usando as estatísticas do conjunto de treinamento)
x_test_features[:, i] = (x_test_features[:, i] - avg) / stddev
print(x_train_features.shape)
print(x_test_features.shape)
Da mesma forma que o exemplo do livro, vamos nos concentrar em identificar apenas os dígitos $1$ e $5$. Vamos juntar os dados com labels 1 e 5, e embaralhar a ordem.
# Antes de mais nada, definir a seed aleatória como uma constante,
# de forma que todos os experimentos obtenham o mesmo resultado
np.random.seed(56789)
x_train_1 = x_train_features[y_train == 1]
x_train_5 = x_train_features[y_train == 5]
y_train_1 = y_train[y_train == 1]
y_train_5 = y_train[y_train == 5]
x_test_1 = x_test_features[y_test == 1]
x_test_5 = x_test_features[y_test == 5]
y_test_1 = y_test[y_test == 1]
y_test_5 = y_test[y_test == 5]
x_train_features_1_5 = np.concatenate([x_train_1, x_train_5], axis = 0)
y_train_1_5 = np.concatenate([y_train_1, y_train_5], axis = 0).astype('float32')
x_test_features_1_5 = np.concatenate([x_test_1, x_test_5], axis = 0)
y_test_1_5 = np.concatenate([y_test_1, y_test_5], axis = 0).astype('float32')
# Considere o digito 1 como a instancia negativa (-1) e o 5 como a positiva (+1)
y_train_1_5[y_train_1_5 == 1] = -1
y_train_1_5[y_train_1_5 == 5] = +1
y_test_1_5[y_test_1_5 == 1] = -1
y_test_1_5[y_test_1_5 == 5] = +1
def shuffle(X, y):
# Com cuidado para embaralhar a entrada e a saída da mesma forma
perm = np.random.permutation(len(X))
return X[perm], y[perm]
x_train_features_1_5, y_train_1_5 = shuffle(x_train_features_1_5, y_train_1_5)
print(x_train_features_1_5.shape, y_train_1_5.shape)
print(x_test_features_1_5.shape, y_test_1_5.shape)
Em seguida, vamos plotar o conjunto de treinamento no plano cartesiano, onde os eixos correspondem aos atributos que acabamos de extrair.
plt.figure(figsize=(20,10))
# Plotamos os 1s em azul...
plt.scatter(x = x_train_1[:,0], y = x_train_1[:,1], c = 'blue', alpha = 0.4)
# ...e os 5s em verde
plt.scatter(x = x_train_5[:,0], y = x_train_5[:,1], c = 'green', alpha = 0.4)
plt.xlabel('Intensidade Média')
plt.ylabel('Simetria')
plt.show()
Como era de se esperar, os $1$s geralmente possuem menor intensidade maior simetria que os $5$s. Além disso, existe mais variação na forma como é possível desenhar os $5$s (note a dispersão dos pontos verdes).
É importante notar que existe uma sobreposição considerável de $1$s e $5$s, indicando a não separabilidade linear dos dados quando representado através destas features.
O primeiro algoritmo avaliado será a versão básica do Perceptron. Iremos rodar o PLA por 1000 iterações, sempre escolhendo o primeiro exemplo erroneamente classificado do dataset.
# Retorna a saida do perceptron definido pelo vetor de pesos w (shape = (d,))
# aplicado aos exemplos na matriz X (shape = (N, d))
def predict(X, w):
# Make sure the data matrix has a bias coordinate
if X.shape[1] != w.shape[0]:
# Add a bias value 1 as the first coordinate of each vector
X = np.concatenate([np.ones((len(X), 1)), X], axis = 1)
return np.sign(np.dot(X, w))
# A função a seguir recebe a matriz de dados X (shape = (N, d)) e vetor de saída
# y (shape = (N,)), e retorna um vetor de pesos de acordo com o PLA
# e, caso return_history = True, uma lista com a quantidade de erros cometidos a cada iteração
def PLA(X, y, w0 = None, max_iterations = 1000, return_history = False):
weights = np.asarray([1, -0.5, 1]) # pesos iniciais
Xe = np.hstack(( np.ones((X.shape[0],1)), X ))
error_history = []
for i in range(1000):
current_predict = predict(X, weights)
sum_errors = np.count_nonzero(current_predict != y)
error_history.append(sum_errors)
first_error_index = (current_predict != y).tolist().index(True)
if y[first_error_index] > 0:
weights += Xe[first_error_index]
else:
weights -= Xe[first_error_index]
return weights, error_history
np.random.seed(1257)
w_pla, errors = PLA(x_train_features_1_5, y_train_1_5, return_history = True)
print(w_pla)
plt.figure(figsize = (12, 8))
plt.plot(errors)
plt.yscale('log')
plt.xlabel('Iteration #')
plt.ylabel('Misclassified Samples')
plt.show()
Observe que a quantidade de erros oscila significativamente ao longo das iterações. Isto se deve a natureza local das atualizações realizadas pelo PLA. Ao corrigir a classificação de um ponto, outros pontos (possivelmente vários) corretos podem passar para o lado errado da fronteira de classificação.
Vamos visualizar a fronteira de classificação obtida junto com os dados.
plt.figure(figsize=(20,10))
# Plotamos os 1s em azul...
plt.scatter(x = x_train_1[:,0], y = x_train_1[:,1], c = 'blue', alpha = 0.4)
# ...e os 5s em verde
plt.scatter(x = x_train_5[:,0], y = x_train_5[:,1], c = 'green', alpha = 0.4)
# A fronteira é uma linha no formato w_pla[0] + w_pla[1]*intensidade + w_pla[2]*simetria = 0
# Obter dois pontos quaisquer na fronteira
p1 = (-2, -(w_pla[0] - 2*w_pla[1])/w_pla[2])
p2 = (2, -(w_pla[0] + 2*w_pla[1])/w_pla[2])
lines = plt.plot([p1[0], p2[0]], [p1[1], p2[1]], '-')
plt.setp(lines, color='r', linewidth=4.0)
plt.xlabel('Intensidade Média')
plt.ylabel('Simetria')
plt.show()
import seaborn as sn
import pandas as pd
y_true = y_train_1_5
y_pred = predict(x_train_features_1_5, w_pla)
true_positives = np.sum((y_pred == +1) * (y_true == +1))
true_negatives = np.sum((y_pred == -1) * (y_true == -1))
false_positives = np.sum((y_pred == +1) * (y_true == -1))
false_negatives = np.sum((y_pred == -1) * (y_true == +1))
confusion = [
[true_positives, false_positives],
[false_negatives, true_negatives]
]
df_cm = pd.DataFrame(confusion, ['$\hat{y} = 5$', '$\hat{y} = 1$'], ['$y = 5$', '$y = 1$'])
plt.figure(figsize = (10,7))
sn.set(font_scale=1.4)
sn.heatmap(df_cm, annot=True, annot_kws={"size": 16}, cmap = 'coolwarm')
Como já era de se esperar pela fronteira de classificação, os $1$s foram classificados mais facilmente do que os $5$s.
Desta vez, vamos usar a versão pocket do Perceptron.
# A função a seguir recebe a matriz de dados X (shape = (N, d)) e vetor de saída
# y (shape = (N,)), e retorna um vetor de pesos de acordo com o Pocket Perceptron
# e, caso return_history = True, uma lista com a quantidade de erros cometidos a cada iteração
def PLA_pocket(X, y, w0 = None, max_iterations = 1000, return_history = False):
weights = np.asarray([1, -0.5, 1]) # pesos iniciais
Xe = np.hstack(( np.ones((X.shape[0],1)), X ))
error_history = [Xe.shape[0]]
print(error_history)
for i in range(1000):
current_predict = predict(X, weights)
first_error_index = (current_predict != y).tolist().index(True)
if y[first_error_index] > 0:
weights += Xe[first_error_index]
else:
weights -= Xe[first_error_index]
sum_errors = np.count_nonzero(current_predict != y)
if sum_errors < error_history[-1]: # Pocket
pocket_weights = weights
error_history.append(sum_errors)
return pocket_weights, error_history
np.random.seed(18375)
w_pocket, errors = PLA_pocket(x_train_features_1_5, y_train_1_5, return_history = True)
print(w_pocket)
plt.figure(figsize = (12, 8))
plt.plot(errors)
plt.yscale('log')
plt.xscale('log')
plt.xlabel('Iteration #')
plt.ylabel('Misclassified Samples')
plt.show()
Como era esperado, a quantidade de erros associados ao vetor de pesos que está "no pocket" jamais aumenta.
Vamos visualizar a nova fronteira de classificação.
plt.figure(figsize=(20,10))
# Plotamos os 1s em azul...
plt.scatter(x = x_train_1[:,0], y = x_train_1[:,1], c = 'blue', alpha = 0.4)
# ...e os 5s em verde
plt.scatter(x = x_train_5[:,0], y = x_train_5[:,1], c = 'green', alpha = 0.4)
# A fronteira é uma linha no formato w_pla[0] + w_pla[1]*intensidade + w_pla[2]*simetria = 0
# Obter dois pontos quaisquer na fronteira
p1 = (-2, -(w_pocket[0] - 2*w_pocket[1])/w_pocket[2])
p2 = (1, -(w_pocket[0] + 1*w_pocket[1])/w_pocket[2])
lines = plt.plot([p1[0], p2[0]], [p1[1], p2[1]], '-')
plt.setp(lines, color='r', linewidth=4.0)
plt.xlabel('Intensidade Média')
plt.ylabel('Simetria')
plt.show()
import seaborn as sn
import pandas as pd
y_true = y_train_1_5
y_pred = predict(x_train_features_1_5, w_pocket)
true_positives = np.sum((y_pred == +1) * (y_true == +1))
true_negatives = np.sum((y_pred == -1) * (y_true == -1))
false_positives = np.sum((y_pred == +1) * (y_true == -1))
false_negatives = np.sum((y_pred == -1) * (y_true == +1))
confusion = [
[true_positives, false_positives],
[false_negatives, true_negatives]
]
df_cm = pd.DataFrame(confusion, ['$\hat{y} = 5$', '$\hat{y} = 1$'], ['$y = 5$', '$y = 1$'])
plt.figure(figsize = (10,7))
sn.set(font_scale=1.4)
sn.heatmap(df_cm, annot=True, annot_kws={"size": 16}, cmap = 'coolwarm')
Desta vez, vamos usar o sinal da regressão linear para classificar os exemplos.
Queremos minimizar o erro quadrado médio:
$$E_{in}(\mathbf{w}) = \frac{1}{N} \sum_{n=1}^N (\mathbf{w}^T \mathbf{x}_n - y_n)^2$$
Para o qual existe uma fórmula fechada (demonstração no livro Learning from Data [2]):
$$\mathbf{w}_{lin} = \text{X}^{\dagger} y$$
onde
$$\text{X}^{\dagger} = (\text{X}^T \text{X})^{-1} \text{X}^T$$
# A função a seguir recebe a matriz de dados X (shape = (N, d)) e vetor de saída
# y (shape = (N,)), e retorna o (único) vetor de pesos que minimiza o erro quadrado médio.
def linear_regression(X, y):
X = np.c_[np.ones(X.shape[0]), X]
w = np.linalg.inv(X.T@X)@X.T@y
return w
w_lin = linear_regression(x_train_features_1_5, y_train_1_5)
print(w_lin)
plt.figure(figsize=(20,10))
# Plotamos os 1s em azul...
plt.scatter(x = x_train_1[:,0], y = x_train_1[:,1], c = 'blue', alpha = 0.4)
# ...e os 5s em verde
plt.scatter(x = x_train_5[:,0], y = x_train_5[:,1], c = 'green', alpha = 0.4)
# A fronteira é uma linha no formato w_lin[0] + w_lin[1]*intensidade + w_lin[2]*simetria = 0
# Obter dois pontos quaisquer na fronteira
p1 = (-2, -(w_lin[0] - 2*w_lin[1])/w_lin[2])
p2 = (1, -(w_lin[0] + 1*w_lin[1])/w_lin[2])
lines = plt.plot([p1[0], p2[0]], [p1[1], p2[1]], '-')
plt.setp(lines, color='r', linewidth=4.0)
plt.xlabel('Intensidade Média')
plt.ylabel('Simetria')
plt.show()
É possível perceber que a fronteira obtida com a regressão linear faz sentido em relação as classes.
import seaborn as sn
import pandas as pd
y_true = y_train_1_5
y_pred = predict(x_train_features_1_5, w_lin)
true_positives = np.sum((y_pred == +1) * (y_true == +1))
true_negatives = np.sum((y_pred == -1) * (y_true == -1))
false_positives = np.sum((y_pred == +1) * (y_true == -1))
false_negatives = np.sum((y_pred == -1) * (y_true == +1))
confusion = [
[true_positives, false_positives],
[false_negatives, true_negatives]
]
df_cm = pd.DataFrame(confusion, ['$\hat{y} = 5$', '$\hat{y} = 1$'], ['$y = 5$', '$y = 1$'])
plt.figure(figsize = (10,7))
sn.set(font_scale=1.4)
sn.heatmap(df_cm, annot=True, annot_kws={"size": 16}, cmap = 'coolwarm')
Neste próximo exercício, vamos tentar aplicar regressão logística. Antes de mais nada, vamos implementar a função de predição, que faz uso da função sigmoid, definida abaixo:
$$\sigma(z) = \frac{1}{1 + e^{-z}}$$
def sigmoid(z):
return 1 / (1 + np.exp(-z))
# Retorna a saida da regressão logística definida pelo vetor de pesos w (shape = (d,))
# aplicado aos exemplos na matriz X (shape = (N, d))
def predict_logistic(X, w):
# Make sure the data matrix has a bias coordinate
if X.shape[1] != w.shape[0]:
# Add a bias value 1 as the first coordinate of each vector
X = np.concatenate([np.ones((len(X), 1)), X], axis = 1)
return sigmoid(np.dot(X, w))
Vamos implementar nossa função de custo. No nosso caso, entropia cruzada.
$$E_{in}(\mathbf{w}) = \frac{1}{N} \sum_{n=1}^{N} \ln(1 + e^{-y_n \mathbf{w}^T \mathbf{x}_n})$$
E o seu gradiente, com relação aos pesos. Sabemos que o gradiente de $E_{in}$ é da seguinte forma:
$$\nabla E_{in}(\mathbf{w}) = - \frac{1}{N}\sum_{n=1}^{N} \frac{y_n \mathbf{x}_n}{1 + e^{y_n \mathbf{w}^T \mathbf{x}_n}}$$
As funções a seguir utilizam várias operações vetoriais em numpy.
def cross_entropy(w, X, y):
return np.mean(np.log(1 + np.exp(-y * np.dot(X, w))))
def cross_entropy_gradient(w, X, y):
N = X.shape[0]
return -np.dot(X.transpose(), y / (1 + np.exp(y * np.dot(X, w)))) / N
Por fim, o treinamento do modelo de regressão logística, usando gradient descent.
# A função a seguir recebe a matriz de dados X (shape = (N, d)) e vetor de saída
# y (shape = (N,)), e retorna um vetor de pesos de acordo com o treinamento da regressão logística
# e, caso return_history = True, uma lista com o valor de cross_entropy a cada iteração
def logistic_regression(X, y, batch_size = None, learning_rate = 1e-2, w0 = None, num_iterations = 1000, return_history = False):
N = X.shape[0]
X = np.c_[np.ones(N), X]
if w0 == None:
w0 = np.random.random(X.shape[1])
if batch_size == None or batch_size > N:
batch_size = N
indices = np.random.permutation(N)
cost_history = np.zeros(num_iterations)
for it in range(num_iterations):
X = X[indices]
y = y[indices]
for i in range(0, N, batch_size):
X_i = X[i:i + batch_size]
y_i = y[i:i + batch_size]
pred = predict_logistic(X_i, w0)
w0 -= (1/N)*learning_rate*(X_i.T@(pred - y_i))
if return_history:
pred = predict_logistic(X_i, w0)
cost_history[it] += (1/batch_size)*(-1*y_i.T@np.log(pred) - (1 - y_i).T@np.log(1 - pred))
if return_history:
return w0, cost_history
return w0
np.random.seed(56789)
w_logistic, loss = logistic_regression(x_train_features_1_5, y_train_1_5, num_iterations = 20000, return_history = True)
print(w_logistic)
plt.figure(figsize = (12, 8))
plt.plot(loss)
plt.xlabel('Iteration #')
plt.ylabel('Cross Entropy')
plt.show()
Como é possível perceber pelo gráfico, o custo aparenta ter estagnado a partir de 10000 iterações. Como a regressão logística produz uma probabilidade, ao invés de visualizar a fronteira de classificação, vamos visualizar a intensidade da saída.
No gráfico abaixo, a intensidade de vermelho corresponde a probabilidade atribuída pela regressão logística para cada ponto. Quanto mais próximo de 1, maior é a confiança da regressão de que o ponto em questão seja um $5$.
plt.figure(figsize=(20,10))
y_pred = predict_logistic(x_train_features_1_5, w_logistic)
plt.scatter(x = x_train_features_1_5[:,0], y = x_train_features_1_5[:,1], c = y_pred, cmap = 'coolwarm')
plt.xlabel('Intensidade Média')
plt.ylabel('Simetria')
plt.show()
Como é possível perceber, pontos nas extremidades do cluster possuem baixa probabilidade (azul escuro, maior confiança de ser um $1$), ou alta (vermelho escuro, maior confiança de ser um $5$), enquanto pontos na intersecção das duas classes possuem uma intensidade intermediária, associada com uma saída próxima de 0.5, idicando que a regressão não diferencia com confiança as classes nesta região.
Para atribuírmos categorias para a saída da regressão logística, precisamos definir um threshold entre 0 e 1, que diferencie as saídas positivas das negativas.
Em situações reais, este threshold deve ser definido de acordo com o custo associado a cada tipo de erro (falso positivo vs falso negativo, lembre-se do exemplo do supermercado e da CIA). Nesta situação, vamos manter o raciocínio simples e assumir um threshold de 0.5.
import seaborn as sn
import pandas as pd
threshold = 0.5
y_true = y_train_1_5
y_pred = predict_logistic(x_train_features_1_5, w_logistic) > threshold
true_positives = np.sum(y_pred * (y_true == +1))
true_negatives = np.sum((1-y_pred) * (y_true == -1))
false_positives = np.sum(y_pred * (y_true == -1))
false_negatives = np.sum((1-y_pred) * (y_true == +1))
confusion = [
[true_positives, false_positives],
[false_negatives, true_negatives]
]
df_cm = pd.DataFrame(confusion, ['$\hat{y} = 5$', '$\hat{y} = 1$'], ['$y = 5$', '$y = 1$'])
plt.figure(figsize = (10,7))
sn.set(font_scale=1.4)
sn.heatmap(df_cm, annot=True, annot_kws={"size": 16}, cmap = 'coolwarm')
[1]: François Chollet and others, Keras, https://keras.io, 2015
[2]: Yaser S Abu-Mostafa, Malik Magdon-Ismail, and Hsuan-Tien Lin, Learning from Data, 2012